ANOVA.yb <- function(df, Col_DV, Col_IV, Col_group = NULL) {
  library(tidyr)
  library(car)
  library(plyr)
  library(dplyr)
  library(agricolae)
  
  if (!is.null(Col_group)) {
    df <- unite(df, "group", all_of(Col_group), sep = "_", remove = F)
    df$group <- as.factor(df$group)
    group <- unique(df$group)
  } else {
    df$group <- as.factor(1)
    group <- unique(df$group)
  }
  
  
  result <- list()
  
  for (i in c(1:length(group))) {
    result[[i]] <- list()
    df1 <- df[df$group == group[i],]
    df1 <- unite(df1, "IV", all_of(Col_IV), sep = "_", remove = F)
    df1$IV <- paste("IV", df1$IV, sep = "_")
    df1$IV <- as.factor(df1$IV)
    IV <- unique(df1$IV)
    
    for (j in c(1:length(Col_DV))) {
      result[[i]][[j]] <- list()
      length2 <- function (x, na.rm=T) {
        if (na.rm) sum(!is.na(x))
        else       length(x)
      }
      datac <- ddply(df1, c(Col_IV, "IV"), .drop=T,
                     .fun = function(xx, col) {
                       c(N    = length2 (xx[[col]], na.rm=T),
                         mean = mean    (xx[[col]], na.rm=T),
                         sd   = sd      (xx[[col]], na.rm=T),
                         min  = min     (xx[[col]], na.rm=T),
                         max  = max     (xx[[col]], na.rm=T)
                       )
                     },
                     Col_DV[j]
      )
      datac$se <- datac$sd / sqrt(datac$N) 
      datac$down <- datac$mean - datac$se
      datac$up <- datac$mean + datac$se
      result[[i]][[j]][[3]] <- datac
      
      model <- paste(Col_DV[j], paste(Col_IV, collapse = "*"), sep = "~")
      normal <- data.frame(group = "", DV = "", IV = "", KS = "", KS_P = "", SW = "", SW_P = "")[-1,]
      for (k in c(1:length(IV))) {
        sample <- df1[df1$IV == IV[k], Col_DV[j]]
        ks <- ks.test(sample, pnorm, mean(sample), sd(sample))
        if (sd(sample) != 0) {
          sw <- shapiro.test(sample)
        } else {
          sw <- list()
          sw[[1]] <- 0
          sw[[2]] <- 0
        }
        normal <- rbind(normal, c(as.character(group[i]), as.character(Col_DV[j]), as.character(IV[k]), 
                                  round(ks[[1]], 4), round(ks[[2]], 4), 
                                  round(sw[[1]], 4), round(sw[[2]], 4)))
      }
      names(normal) <- c("group", "DV", "IV", "KS", "KS_P", "SW", "SW_P")
      result[[i]][[j]][[1]] <- as.data.frame(leveneTest(eval(parse(text = model)), data = df1, center = mean))
      result[[i]][[j]][[1]] <- cbind(rep(group[i], length(result[[i]][[j]][[1]][,1])),
                                     rep(Col_DV[j], length(result[[i]][[j]][[1]][,1])),
                                     result[[i]][[j]][[1]])
      names(result[[i]][[j]][[1]])[1:2] <- c("group", "DV")
      result[[i]][[j]][[2]] <- as.data.frame(summary(aov(eval(parse(text = model)), data = df1))[[1]])
      result[[i]][[j]][[2]] <- cbind(rep(group[i], length(result[[i]][[j]][[2]][,1])),
                                     rep(Col_DV[j], length(result[[i]][[j]][[2]][,1])),
                                     row.names(result[[i]][[j]][[2]]),
                                     result[[i]][[j]][[2]])
      names(result[[i]][[j]][[2]])[1:3] <- c("group", "DV", "source")
      SNK <- SNK.test(aov(eval(parse(text = paste(Col_DV[j], "IV", sep = "~"))), data = df1), "IV")$group
      SNK <- data.frame(IV = rownames(SNK), SNK = SNK$groups)
      LSD <- LSD.test(aov(eval(parse(text = paste(Col_DV[j], "IV", sep = "~"))), data = df1), "IV", p.adj="none")$group
      LSD <- data.frame(IV = rownames(LSD), LSD = LSD$groups)
      duncan <- duncan.test(aov(eval(parse(text = paste(Col_DV[j], "IV", sep = "~"))), data = df1), "IV")$group
      duncan <- data.frame(IV = rownames(duncan), duncan = duncan$groups)
      HSD <- HSD.test(aov(eval(parse(text = paste(Col_DV[j], "IV", sep = "~"))), data = df1), "IV")$group
      HSD <- data.frame(IV = rownames(HSD), HSD = HSD$groups)
      mc <- merge(merge(SNK, LSD, by = "IV"), merge(duncan, HSD, by = "IV"), by = "IV")
      result[[i]][[j]][[3]] <- merge(result[[i]][[j]][[3]], merge(normal, mc, by = "IV"), by = "IV")
      result[[i]][[j]][[3]] <- result[[i]][[j]][[3]][, c("group", "DV", Col_IV, "N", "mean", "sd", "min", "max", "se", "down", "up", "KS", "KS_P", "SW", "SW_P", "SNK", "LSD", "duncan" ,"HSD")]
    }
  }
  return(result)
} #方差分析；输入(列名)：df-数据表(数据框)、DV-因变量(向量)、VI-自变量/固定因子(向量)、group-多重比较分组(向量)
#输入：result[[group]][[DV]][[1-levene(n>50-KS); 2-anova; 3-summary & normality_test & multiple_comparisons]]


#示例
data <- data.frame(IV_A = as.factor(rep(c(rep(1,20),rep(2,20),rep(3,20),rep(4,20)),2)),
                   IV_B = as.factor(rep(rep(c(rep(1,5),rep(2,5),rep(3,5),rep(4,5)),4),2)),
                   re = rep(rep(c(1:5),16)),
                   group = as.factor(c(rep(1,80),rep(2,80))),
                   DV_1 = runif(160),
                   DV_2 = rnorm(160))
                           #创建示例数据框，IV_A&IV_B(自变量A、B, 需为因子factor), re(重复), group(分组), DV_1&DV_2(因变量1、2)
data                       #查看示例数据框

result <- ANOVA.yb(data, c("DV_1", "DV_2"), c("IV_A", "IV_B"), c("group"))
                           #分析
result[[1]][[1]][[1]]      #分组1，DV_1，的levene
result[[1]][[2]][[2]]      #分组1，DV_2，的方差表
result[[1]][[1]][[3]]      #分组1，DV_1，的统计描述、正态性检验&多重比较
result[[2]][[2]][[3]]      #分组2，DV_2，的统计描述、正态性检验&多重比较